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Abstract. In this work the Gardner problem of inferring interactions and fields for an 
Ising neural network from given patterns under a local stability hypothesis is addressed 
under a dual perspective. By means of duality arguments an integer linear system is 
defined whose solution space is the dual of the Gardner space and whose solutions 
represent mutually unstable patterns. We propose and discuss Monte Carlo methods 
in order to find and remove unstable patterns and uniformly sample the space of 
interactions thereafter. We illustrate the problem on a set of real data and perform 
ensemble calculation that shows how the emergence of phase dominated by unstable 
patterns can be triggered in a non-linear discontinuous way. 


PACS numbers: 84.35+i, 05.10Ln, 02.40Ft 

1. Introduction 

The use of concepts, methods and models from statistical mechanics, in particular of 
disordered systems and spin glasses, for the analysis of neural networks P] and associative 
memories has a long standing tradition and has given many fruitful insights [2j [,3J |2J [5j |6] . 
Recent years witnessed a renewed surge of interest in modeling real biological neural 
networks with statistical mechanics models since it has been shown that experimental 
observations on the statistics of spiking patters can be reproduced within maximum- 
entropy Ising modelsjT]. The computationally difficult inverse problem of inferring 
interactions and fields from patterns statistics has been addressed in several works([8] 
and references therein), revealing many subtleties concerning for instance the assessment 
of the extrapolated criticality of the inferred model[9j [TOj. In general the problem of 
finding the right parameters of the model in order to store given patterns and structure 
them in an associative memory is an interesting problem of neural network learning. 
In a different and simpler setting with respect to the aforementioned inverse problem, 
interactions and fields could be inferred under a local stability hypothesis of the patterns, 
that is the Gardner problem we will afford herejllj. Consider the patterns = ±1, 
where /i — 1... M is the pattern index and i — 1... N is the neuron index. The pattern 
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H is stable if interactions Jy(we consider a symmetric model Jij = Jji without self 
interactions J lt = 0) and fields h t are such that the system of inequalities 

+ hi)>0 Mi (1) 

jV* 

has solutions (the Hopheld choice J\j = CniCfij K = gives a straightforward 
solution). All the patterns are mutually stable if the inequalities are verified Mp (in 
this case the Hopheld choice is not guaranteed to give a solution). This is a linear 
system of inequalities for the J t j and hf 

'y ' T ^ 0 V/i, i. (2) 

j¥=i 

The problem of finding interactions and fields and characterize statistically the space 
for given ensembles of patterns has been addressed in several works rara ra • 

In this work the Gardner problem is afforded with a dual perspective. From duality 
arguments we define an integer linear system whose solutions give a certificate of the 
infeasibility of the Gardner problem and represent sets of mutually unstable patterns. 
In the following we discuss Montecarlo methods in order to find such sets and sample 
the space of interactions thereafter, we analyze with them a simple instance of biological 
data and perform ensemble calculations. We finally conclude summarizing the results 
and drawing some future perspectives on this problem. 

2. Results 

2.1. The dual Gardner problem 

We have seen that a local stability hypothesis for the patterns imply that the interactions 
and fields of an Ising neural network should verify a system of linear inequalities. The 
existence of solutions for a system of linear inequalities is ruled by duality theorems 
that state that the system has solutions if and only if a suitably defined dual system 
has no solutions jT5]. The space of interactions and holds, system (2), is defined by NM 
inequalities for N(N + l)/2 variables. According to the Gordan theorem the dual of 
system (2) is 

y ' 0 Vi < j (3) 

u 

y ' 0 Vi 

k^i > 0 

that is a system of N(N + l)/2 linear equations for the NM non negative variables 
. In the appendix we report a self consistent demonstration of the theorem based 
on the Fourier-Motzkin-Chernikova elimination method. It should be noticed that the 
matrices defining the constraints are integers and thus any solution can be represented 
in terms of integer solutions, and in essence we can consider system (3) as an integer 
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linear system. Further we note that the equations do not mix variables k^i with different 
% and j. In the following section we will discuss numerical methods to find solutions for 
both systems, the space of interactions and fields and its dual. 

2.2. Resolution methods 

A system of linear inequalities like system ( 2 ) can be solved by relaxational algorithms. 
Starting from a prior (Jk, h°), we calculate the least unsatisfied constraint from (2) and 
update the parameters along the orthogonal direction, defining the series 

(p t , i t ) = argmin^ J^i^j + h%i (4) 

44 - J tj + iZutittiHj K x - K + t Uu ( 5 ) 

that will converge in polynomial time to a solution, provided its existence. The length 

7 of the step can be constant (MinOver[l3]) or proportional to the amount by which 

the constraint is violated (Motzkin algorithm[TBj). In absence of solutions the learning 
through relaxation does not converge, and according to the Gordan theorem there should 
be solutions to the dual system (3). On the other hand, the solutions of the dual system 
(3) can be seen as the ground state with zero energy of the Hamiltonian over integer 
variables obtained summing over the constraints square: 

H = v £(£ *i»u 2 (6) 

i<j [i i 

In this way it is possible to use computational methods from statistical mechanics to 
find solutions. For instance, upon introducing a fictitious inverse temperature [5 we 
can sample equilibrium configurations from the Boltzmann distribution P oc e~^ H 
with a Metropolis-Hastings Montecarlo Markov chain and decrease the temperature 
(simulated annealing) till convergence to a ground state configuration (eventually the 
trivial solution). In order to reduce the dimension of the space we can couple the 
methodspT], be. we can search for solutions of system (3) only among the variables 
that correspond to most unsatisfied constraints for a previous search of solutions to 
system (2). In this way we will have either a solution to the Gardner problem or a set of 
mutually unstable patterns that proves that the Gardner problem is unfeasible for such 
patterns. In the latter case we could reject and remove one of the pattern in the set 
and start again. The choice of the pattern to remove in the set should depend on the 
specific context of the problem afforded. If, on the other hand, the algorithm converge 
to a solution of the Gardner problem, a complete characterization of the solution space 
of system ( 2 ), that is geometrically a polyhedral cone, would require the calculation of a 
so called Hilbert basisp2J, i.e. the complete set of solutions that cannot be represented 
as convex sums of other solutions. This is a difficult computational problem and an 
alternative is to characterize the space spanned by solutions of system ( 2 ) with statistical 
methods. This would require to fix a scale for the space and a simple choice consists 
in bounding interactions and fields in way that depends on the specific context of the 
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problem afforded, for instance Rj G [—C,C\ hi G [—(7,(7]. In this way the space of 
interactions form a convex polytope. An uniform sampling of convex space in high 
dimensions can be performed in feasible times by means of Monte Carlo Markov chains, 
the Hit and Run method being the faster so far1 19] [20]. Given a /4-dimensional convex 
set P, from which one wants to sample from, and a point inside x t G P, the Hit and 
Run algorithm is defined as follows: 

(i) Choose a uniformly distributed direction 9 t , that is, a point extracted from the 
uniform distribution on the 14-dimensional unit sphere. This can be done with the 
Marsaglia method, i.e. by generating D independent gaussian random variables 9 l t 
with zero mean and unit variance, and then normalizing the vector to unit length; 

(ii) Extract A* uniformly from the interval [A min , A max ], where A min (A max ) is the 
minimum (maximum) value of A such that Xt + A 9 t G P; 

(iii) Compute the point x t+ i — x t + A *9 t , increment t by one and start again. 

The mixing time r of the hit and run scales like [2D] 

r ~ 0{D 2 R 2 /r 2 ) (7) 

where D is the dimension of the polytope, P, r are the radii of respectively the minimum 
inscribing and the maximum inscribed balls. The factor R/r can be large and leading 
to ill-conditioning but it can be reduced to y/D for centrally symmetric polytopes 
and to D in general, by an affine transformation defined by the so-called Loewner- 
Jolin Ellipsoidal]. i.e. the ellipsoid of maximal volume contained in the polytope. 
Unfortunately this ellipsoid cannot be found in polynomial time, but it has been shown 
by L. Lovazs that a weaker form of the Loewner-John ellipsoid, with a factor of P 3//2 , 
can be found in polynomial time [22]. 

2.3. Analysis of real data: an illustrative example 

In this section we will analyze a subset of experimental data on the neurons spiking 
patterns in a salamander retina, in particular we refer to the snapshot of T = 2 s 
reported in Fig 1 of[23]. There are M — 31 distinct patterns for N = 15 effectively 
spiking neurons that we report in the appendix. The relaxational algorithm applied 
to the space of synaptic strengths for this network does not converge indicating the 
presence of conflicting patterns that have been identified with a Montecarlo simulated 
annealing. For sake of simplicity, we have chosen to remove one of the pattern picked 
uniformly at random in the retrieved set. We report in the appendix the set of conflicting 
patterns retrieved in one run of the algorithm and the patterns we have removed. Once 
we have removed a sufficient number of conflicting patterns such that the system (2) 
defining the space of synaptic strength is feasible we can sample interactions and fields. 
We have chosen to bound G [—1000,1000] hi G [—1000,1000] and we have sampled 
uniformly the space by means of the Hit-and-run Markov chain after ellipsoidal rounding 
preprocessing]!?!]. In Fig 1 we show the length of the diameters of the ellipsoid that 
rounds the polytope (the dimension is D = N(N + l)/2 = 120) for decreasing order, 
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Figure 1 . diameters length of the ellipsoid that rounds the polytope representing the 
space of neural network models storing a given set of stable patterns (see appendix). 

where we can appreciate the heterogeneous structure of the space (the ratio between 
the largest and smallest diameter is of order 10 4 ). In Fig 2 we show the averages of the 



Figure 2. Average interaction strengths and fields retrieved from uniform sampling of 
the space of neural network models corresponding to the stable set of patterns reported 
in the appendix. 


interaction strengths (off-diagonal) and fields (diagonal) thus retrieved. Fields are on 
average negative for all neurons apart from the fourth that has positive (ferromagnetic) 
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interactions with the other neurons (whose interactions are on average weakly anti 
ferromagnetic) 


2 . 4 . Ensemble calculations in a simple setting 


In order to get a picture of the dependence on N and M of the volume of the solution 
space of (3), the problem is that, at odds with linear algebra, for integer linear systems 
there are no straightforward results relating it for instance to the kernel of the constraints 
matrix. In order to get such insights we will perform here a statistical analysis of the 
solution space of (3) by means of ensemble calculations in a simple setting. We will 
consider for simplicity an integer linear system for the variables k^ = 0 , 1 : 

Y k ^i = 0 ( 8 ) 

u 

with — M , % — 1... N . We want to enumerate solutions for random systems with 
i = ±1 with equal probability (a maximum entropy case). They are the ground states 
with zero energy of the system with hamiltonian 

H = ^E<E«*<) 2 ' w 

i [i 

We would like to evaluate log Z, where Z = )T) fc =0 , e~ /3H and the bar stands for the 
average over disorder. We will calculate Z n , and perform the limit 

log Z = lim — log Z n . (10) 

n->0 n 

Then we will send f3 —> 00 , and N, M —> 00 with a = finite. We have (a — 1 ... n is 
the replica index) 


Z n 


eii e -/3/;v(E M b-bu) 2 = 

JO* 


=zn 


3 P Ea biaki/a( Ej j 


( 11 ) 


where we have simply developed the square and exchanged index in the sums and 
products. Now we observe that for N large 


N 


i&j 


1 if fi — v , 

->■ x^„ if 


( 12 ) 


where are independent gaussian random variables centered in 0 with std . Then 
we have 


tr J ti -znx 

=En g/3 2 /V(E a kliakva) 2 n g ft Ea k/iakjja __ 


N 2 g- P Ea = 




En g/3 2 W((E,j — E m n e ^ 




(13) 
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where in the last line we have used k= k^a, however the term for the correction p ^ v 
is not extensive and we will drop it in the following. We use known property of Gaussian 
integrals: 


e i(^E^GW 2 = / d,Q ab 


e ~N^Db+pQ ab 




(14) 


we have then 


Z n = 






a, b 


e~ N ^ 


y/4n/N 


lit 

k a a,6 


,/3k a Qabk b J~~J" e ~pk a ^ M _ 


= / dQe MF W) 


F(g) = _ 4 ^S g ^ + logZr(Q) 


a, b 


Z r = J2 e ^ Hr H r = -J2kaQabk b + J2 k - ( 15 ) 


k a ab a 

We see that saddle point equations lead to the following interpretation of the Q ab : 
OF 


dQab 


= 0 


Off = 2 af)(k a k b ) 


(16) 


where the brackets stand for averages among interacting replicas with hamiltonian H r . 
A natural ansatz is to assume replica-symmetric matrices 


Qab — 


q 0 if a = b , 


qi if a ^ b . 
and we have the expression 


(17) 


n 


F(qo,qi) = ln i'/rl + (n - i )ql) + log z r 


Z r — — 


d,Xe~^~ (1 + e vm7\-d(i+qi-q 0 )^ 


(18) 


where we have been using again known properties of gaussian integrals. The free energy 
/ in the limit n —> 0 


— Pf — lim F/n = 

17—>0 

= - — (<?0 - dl) + —L [ dXe~^ log(l + e^-W+51-®)) 
4 a v 2tt J 

If we analytically continue the SP equations for n —> 0 we obtain: 

i r „ a 2 i 


(19) 


q 0 = 2a/3- 


Qi = Qo~ 2ol\ 


A 

dXe~^ 


l _|_ e -\/2pqi\+P(l+qi-qo) 

P 1 f .. A 


2 gi 


dAe 2 


X -)- g — p2/3giA+^(l+9i— go) 


( 20 ) 
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In the limit f3 —> oo we look at the solution 


lim q 0 /j3 — lim q\//3 — x 

/3—>oo j3—>oo 


x = 2 a 


f = 


y/2-n 
1 


d\e 


' 1 
\/2x 
f* OO 




> i 
y/2x 


d\e 2 (-\/2xA — 1) 


( 21 ) 

( 22 ) 


We remind that the original problem is for NM variables and N(N — l)/2 equations 
and thus we should consider a/2 as the correct control parameter, further upon looking 
at the interpretation of the saddle point equations, we will consider r = 2 ax as an order 
parameter in (0,1) Apart from the trivial solution r — 0 we depict in the following figure 
the curves r(a ) and f(a). 




Figure 3. Order parameter r and (minus) free energy / as a function of a = M/N 
from replica symmetric ensemble calculations. 


3. Conclusions 

In this work we have defined the dual of the space of interactions of neural network 
models under a stability hypothesis for the patterns. By means of the Gordan theorem 
we derived an integer linear system whose solutions represent mutually unstable patterns 
whose feasibility implies the infeasibility of the Gardner problem. We have discussed 
Monte Carlo methods for detecting and removing such instabilities and finally perform 
an uniform sampling of the space of interactions and fields. We have applied them 
for illustrative purposes to the analysis of a simple data set of spiking patterns from 
biological neurons. We have performed replica-symmetric ensemble calculations in a 
simple setting showing that the emergence of a phase dominated by unstable patterns 
can be triggered in a discontinuous way. Apart from an analysis of more complex 
ensembles it would be important to perform stability analysis, finite temperature 
calculations and for for general integer variables /q = 0,1,..., m. Recent analysis 
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seems to show however that for sparse random integer linear system the behavior for 
m > 1 is qualitatively similar to that for rri — 1 and that the replica symmetric solution 
is stable [25]. Regarding applications, for instance in analyzing real data of spiking 
neurons, we stress that our choice to remove a pattern uniformly at random among the 
ones retrieved in unstable sets is motivated by simplicity and for illustrative purposes. In 
order to implement different choices it would be worth to have an overall picture of the 
space of mutually unstable patterns for single instances, beyond our Montecarlo method 
that retrieves just one solution. Probably this problem could gain fruitful insights from 
message passing and cavity methods, that have been recently applied to integer linear 
systems but in a different context [22] [27] [28]. 


Appendix 


Demonstration of the Gordan theorem 

Consider the system of homogeneous linear inequalities for the variables x M defined by 
the matrix A- 




/ J Ai^x^ 0 


Vi. 


(23) 


. We associate it to the following dual system for the non-negative variables yt 
^ 0 V fi 

i 

Vi> 0, y ^ 0 

We want to demonstrate the Gordan theorem: 


(24) 


One and only one of the systems (23) and (24) have solution. 
We first demonstrate that 


(24) has solution (23) has no solution. 


Suppose there is a solution y* of (24). Take any vector x, multiply it component by 


component with the equations of the system (24) and sum over them. Exchanging the 
indeces p and i in the sums and given the positivity of y*, it is straightforward to 


conclude that no vector x can satisfy the system (23). 
Now we demonstrate that 


(23) has no solution => (24) has solution. 


The demonstration is given by induction in the number M of unknowns x^. The 


statement is true for M — 1. Infact, for one unknown, the system (23) is inconsistent if 


and only if there is at least one couple i, j such that AiA'i = ~~ 1/c < 0; and in this case 


y, = 1, yj = c and yi = 0 VI ^ i,j is a solution of (24). Let’s consider a system of the 


type (23) with M unknowns and suppose it is inconsistent. We will prove that the dual 
system has solutions supposing the theorem true for systems with M — 1 unknowns. We 
have Vi Y^l=\ > - Am^m, if Am 7^ 0, we can define A ifl = — A^/Am, and we 
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have: 

M-1 

^ ' Ai^X^ Pi > iCM Vi • '4 v A/ “C 0 

M 

M—l 

^ ' A MV Qj ^ .r: a/ Vj . AjM x 0 

M 

M-l 

Y J K x ^ = R l > 0 VZ: Am = 0 . ( 25 ) 

M 

Writing the system in this form, we can pass to the following system in M — 1 unknows: 


Pi > Qj 

Ri > 0 


Vi, j : Am < 0 
VZ : Am = 0. 


Am > 0 


(26) 


Now, this system is also inconsistent. 

Suppose infact there is a solution x* = (x*), /a = 1... M — 1. 

We could add to it any x* M such that maXjQj(pt*) < x* M < miniPiipt*) and we will have 
a solution for the original system as well, against the hypothesis. By induction, the 
theorem is true for systems with M — l unknowns. Then, referring to (26), there are 
Vij > 0, yi > 0 with at least one positive, such that 'YhijVij (Am — Am) + Xu ViA-in - 
0 V/r. 

From this we have finally a solution for the system (24): 

Vi 


Vi / J Vij /A-iM 

3 

Uj ~ 'y ] Vij/AjM 


Am < 0 


Vi : Am > 0 


yi 


VZ : Aim = 0, 


(27) 


and the theorem is proven. Regarding the application to the Gardner problem it is 
useful to remind that, essentially if we group index r = (p,i), s = (i, j) the matrix is 


Am4)6j) 


A*(0) with (without) fields if i = j , 

if *7A ■ 


(28) 
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We report here the patterns analyzed in sec 2.3 taken from | 23 j . i.e. the 0 — 1 matrix 
where = 1 (0) means that the neuron i is firing (not) in the pattern p. 


NEURONS 



0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 
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0 

0 
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1 

0 
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0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

p 

A 

T 

T 

TP 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 
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0 

0 

0 

1 

0 

0 

0 
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0 
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0 
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0 

0 
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0 

0 

0 

0 

0 

0 

1 
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0 
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0 
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1 

0 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

0 


We have been using the usual transformation £,; /t = 2 — 1 . We report in the following 
table 1 the unstable patterns solutions of system ( 3 ) retrieved and removed in one run 
of the algorithm described in sec 2 . 2 . The information given in the table should be 
interpreted in the following way: if patterns 29 and 6 are conflicting for neuron 4 , 
this means that system ( 3 ) has a solution k 4^9 = = 1 and = 0 for the other 

components. 
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Conflicting patterns 

Neuron 

Pattern removed 

29,6 

4 

29 

30,31 

7 

31 

2,5 

11 

5 

9,10 

13 

9 

1 , 10 , 11,14 

7 

14 

16,26 

13 

26 

20,23 

1 

23 

7,28 

3 

28 

8,16 

3 

8 

1 , 6 , 7,25 

12 

6 

4,19 

4 

19 

1 , 2 , 20,24 

11 

24 

21,22 

11 

21 

3,20 

1 

3 


Table 1. First Column: Set of conflicting patterns. Second Column: Reference 
neuron. Third column: Removed pattern. Obtained by coupling Relaxation and 
Montecarlo simulated annealing applied to the data represented in matrix (20). 
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